library(lmtest)
library(sandwich)
library(date)
library(doBy)



load(file="data/sb_2.Rdata")##load full data with block groups coded
dim(sb)

sb<-subset(sb, sb$suspected.crime=="cpw")
dim(sb)

sb<-sb[order(sb$dayno),]

sb$wpct<-as.numeric(as.character(sb$white))/as.numeric(as.character(sb$total))
summary(sb$wpct)

#load most recent version of median
load("data/med.white.Rdata")##load median % white bg (constructed in data management file)
cut<-med.white


memo<-mdy.date(year=2013, month=3, day=5)

sb$b<-NA
sb$b[sb$wpct<cut]<-0
sb$b[sb$wpct>=cut]<-1 ##these are very white places
table(sb$b)

sb<-sb[order(sb$dayno),]





#CONSTRUCT LAGGED DV
means<-summaryBy(weapon.hit~dayno, data=sb, FUN=mean, na.rm=T)##get daily hit rates
means$lag.hit.new<-NA
means<-means[order(means$dayno, decreasing=F ),]
head(means)

for(i in 2:length(means$dayno)){
	
	means$lag.hit.new[i]<-means$weapon.hit.mean[i-1]
	
}


head(means)


dim(sb)
sb<-merge(sb, means, by="dayno",all.x=T)
head(sb)
dim(sb)

sb<-sb[order(sb$dayno),]



sb$weekday<-weekdays(as.Date(sb$date))


sb$period<-NA
sb$period[sb$dayno>=memo]<-"post-memo"
sb$period[sb$dayno<memo]<-"pre-memo"
sb$period<-factor(sb$period, levels=c("pre-memo","post-memo"))
table(sb$period)

sb$distance<-sb$dayno-memo

sw2<-sb

sw2$int<-sw2$b


##TABLE 1
## Estimate discontinuity at moment of intervention using all weapon stops, 2008-2015

se1<-rep(NA, 8)
names.models<-c("diff.means","diff.means+controls","linear","linear+controls","squared","squared+controls","cubic","cubic+controls")
names(se1)<-names.models

##diff in means
m1<-lm(weapon.hit~period+int+period*int, data=sw2)
se1[1]<-summary(m1)$coefficients["periodpost-memo:int","Std. Error"]

m1a<-lm(weapon.hit~period+int+period*int+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=sw2)
se1[2]<-summary(m1a)$coefficients["periodpost-memo:int","Std. Error"]

##local linear
m2<-lm(weapon.hit~period+distance+period*distance + int+period*int+distance*int+period*distance*int, data=sw2)
se1[3]<-summary(m2)$coefficients["periodpost-memo:int","Std. Error"]


m2a<-lm(weapon.hit~period+distance+period*distance + int+period*int+distance*int+period*distance*int+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=sw2)
se1[4]<-summary(m2a)$coefficients["periodpost-memo:int","Std. Error"]

##second order poly
m3<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2), data=sw2)
se1[5]<-summary(m3)$coefficients["periodpost-memo:int","Std. Error"]

m3a<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=sw2)
se1[6]<-summary(m3a)$coefficients["periodpost-memo:int","Std. Error"]


##cubic
m4<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2)+int*I(distance^3)+int*period*I(distance^3), data=sw2)
se1[7]<-summary(m4)$coefficients["periodpost-memo:int","Std. Error"]

m4a<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3)+int+period*int+distance*int+I(distance^2)*int+period*int*distance+period*int*I(distance^2)+int*I(distance^3)+int*period*I(distance^3)+factor(year)+factor(month)+factor(weekday)+lag.hit.new, data=sw2)
se1[8]<-summary(m4a)$coefficients["periodpost-memo:int","Std. Error"]


models2<-list(m1, m1a, m2, m2a, m3, m3a, m4, m4a)

print("models done")

save(se1,file= "output/se1_white_bg.Rdata")##save conventional SEs
save(models2,file= "output/std_models_white_bg.Rdata")##save conventional SEs

vars<-list(NA)
vc<-list(NA)

##estimate HAC SEs


print("computing HAC SEs")

for(i in 1:length(models2)){
vci<-vcovHAC(models2[[i]])
vc[[i]]<-vci
save(vc, file="output/white_bg_HAC_vc.Rdata")
print(i)
}
